function [t_Ws_L,S_Ws_L,...
          t_Ws_R,S_Ws_R,...
          t_Wu_R,S_Wu_R,...
          t_Wu_L,S_Wu_L] = PropagateManifold_Colinear(~,varargin)

%%

numvarargs = length(varargin);
% Number of optional arguments

if mod(numvarargs, 2) == 1
    error('Error: Please check name-value pair arguments');
end

name_arguments = varargin(1:2:end);
value_arguments = varargin(2:2:end);
% Initialize name-value arguments

for ii = 1:length(name_arguments)
    switch lower(name_arguments{ii})
        case 'poinitialstate'
            Xo = value_arguments{ii};
        case 'poperiod'
            tau = value_arguments{ii};
        case 'massratio'
            u = value_arguments{ii};
        case 'manifoldtrajectoryamount'
            numpoints = value_arguments{ii};
        case 'integrationtime'
            integrationTime = value_arguments{ii};
        case 'system'
            if value_arguments{ii} == 'SE'
                epsilon=5*10^(-6);
            elseif value_arguments{ii} == 'EM'
                epsilon=5*10^(-5);
            end

    end
end


%% Initial Guess
% northern

% tau = [1.39253390838970];
%[] Orbital period initial guess.
T = tau;
%[] Single orbital period

to = linspace(0,T,numpoints);
%[]

options=odeset('RelTol',1e-15,'AbsTol',1e-21);
%[] Increase the Ode 45 tol

[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S';
t = t';
%[] Integrate states.

Monodromy = State2STM(t,S,u);
%[] Monodromy matrix

[eig_vec,eig_val] = eigs(Monodromy);
%[] Find eigen values and vectors

eig_val = diag(eig_val);
%[] Get just the value.

[~,ind_stable] = min(eig_val);
[~,ind_unstable] = max(eig_val);
%[] Determine stable and unstable index

nu_stable = eig_vec(:,ind_stable);
nu_unstable = eig_vec(:,ind_unstable);
%[] Stable and unstable direction perturbation

InitialConditions_stable_outside = zeros(6,length(t));
InitialConditions_stable_inside = zeros(6,length(t));
InitialConditions_unstable_outside = zeros(6,length(t));
InitialConditions_unstable_inside = zeros(6,length(t));

%[] Pre allocate initial conditions matrix

% Loop for creating initial conditions for manifold trajectories.
parfor ii = 1:(length(t)-1)
    
    t_STM = [t(1),t(ii+1)];
    %[] Determine the state transition matrix at desired time T
    
    Pert_unstable = State2STM(t_STM,Xo,u)*nu_unstable;
    Pert_stable = State2STM(t_STM,Xo,u)*nu_stable;
    %[] calculate the perturbance needed for unstable and stable manifold
    %trajectories
    
    Norm_Pert_unstable = Pert_unstable / norm(Pert_unstable);
    Norm_Pert_stable = Pert_stable / norm(Pert_stable);
    %[] Normalize the perturbing force.
    
    if Norm_Pert_stable(1)<0 % If x perturbation direction is negative
        
        InitialConditions_stable_outside(:,ii)      = S(:,ii) + epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) - epsilon * Norm_Pert_stable;
%         disp('a')
    else
        InitialConditions_stable_outside(:,ii)      = S(:,ii) - epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) + epsilon * Norm_Pert_stable;
%         disp('b')
    end
    
    if Norm_Pert_unstable(1)<0 % If x perturbation direction is negative
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) + epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) - epsilon * Norm_Pert_unstable;
%         disp('a')
    else
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) - epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) + epsilon * Norm_Pert_unstable;
%         disp('b')
    end
    
    %[] Determine the initial states for stable manifold's initial states.
    
end
T_manifold_stable_outside = integrationTime:-0.0001:0; % Integrating backward in time
%[] Integration time for stable manifold coming from outside

T_manifold_unstable_inside = 0:0.00025:integrationTime;% Integrating forward in time
% options_SecondBodyCrossing=odeset('RelTol',1e-6,'AbsTol',1e-9,'Events',@(t,S)HaltProp(t,S,u));

%Loop for propagating the calculated initial conditions
parfor ii = 1:length(InitialConditions_stable_outside)-1
   
    [t_Ws_L{ii},S_Ws_L{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_stable_outside,...
        InitialConditions_stable_outside(:,ii),options);
    S_Ws_L{ii} = S_Ws_L{ii}';

%     if sum(S_stable_outside_SM{ii}(1,:)>(1-u))>0
%         S_stable_outside_SM{ii} = nan;
%     end
    %[] Integrate the stable manifold trajectory coming in from outside.
%     
    [t_Ws_R{ii},S_Ws_R{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_stable_outside,...
        InitialConditions_stable_inside(:,ii),options);
    S_Ws_R{ii} = S_Ws_R{ii}';
    %[] Integrate the stable manifold trajectory coming in from inside.
%     
    [t_Wu_R{ii},S_Wu_R{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_unstable_inside,...
        InitialConditions_unstable_inside(:,ii),options);
    S_Wu_R{ii} = S_Wu_R{ii}';
    %[] Integrate the stable manifold trajectory coming in from inside.
    
    [t_Wu_L{ii},S_Wu_L{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_unstable_inside,...
        InitialConditions_unstable_outside(:,ii),options);
    S_Wu_L{ii} = S_Wu_L{ii}';
    
%     if sum(S_unstable_outside_SM{ii}(1,:)>(1-u))>0
%         S_stable_outside_SM{ii} = nan;
%     end
    %[] Integrate the stable manifold trajectory coming in from inside.
    
end

end
